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Fine structure generation in double-diffusive system 
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Double-diffusive convection in a horizontally infinite layer of a unit height in a large Rayleigh 
numbers limit is considered. From linear stability analysis it is shown, that the convection tends 
to have a form of travelling tall thin rolls with width about 30 times less than height. Amplitude 
equations of ABC type for vertical variations of amplitude of these rolls and mean values of diffu- 
sive components are derived. As a result of its numerical simulation it is shown, that for a wide 
variety of parameters considered ABC system have solutions, known as diffusive chaos, which can 
be useful for explanation of fine structure generation in some important oceanographical systems 
like thermohaline staircases. 
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I. INTRODUCTION 



Double-diffusive or thermohaline convection plays im- 
portant role in a heat-mass transfer processes in the 
ocean p|. It also essentially influences on different small 
scale processes, like formation of vertical temperature 
and salinity fine structure. Such phenomena are not well 
understood at present day. There are only few works, de- 
voted to analytical models of fine structure generation in 
the sea [j, tJ |4| . But no one of them essentially considers 
the role of double diffusion in such processes. As an ex- 
ception one can mention work 5] , where double-diffusive 
step-like fine structure is simulated by numerical Monte 
Carlo methods. The purpose of this article is to develop 
one more mathematical model of two dimensional double- 
diffusive convection in a horizontally infinite layer, based 
on system of amplitude equations, describing formation 
of vertical fine structure, which in some aspects resem- 
bles actual experimental and observational data. The 
main idea of this work consists in combining the result, 
that in the limit of large Rayleigh numbers convective 
cells tend to be narrow and tall, with constructing of 
ABC-system of amplitude equations @ with respect to 
vertical coordinate in the case of such cells. So this ar- 
ticle includes two main sections apart from this one. In 
the section ITT1 from linear stability problem we determine 
sizes of the most prominent cells in large Rayleigh num- 
bers limit. In the section llTTl in weakly nonlinear approxi- 
mation by multi-scale decomposition technique we derive 
the ABC-system. Its numerical simulation gives vertical 
fine structure via so called diffusive chaos solutions. 

The initial equations describe two-dimensional thermo- 
haline convection in a liquid layer of thickness h, bounded 
by two infinite plane horizontal boundaries. The liquid 
moves in a vertical plane and the motion is described by 
the stream function ip(t, x, z). The horizontal x and ver- 



tical z space variables are used; the time is denoted by 
t. It is assumed, that there are no distributed sources 
of heat and salt, and on the upper and lower bound- 
aries of the area these quantities have constant values. 
Hence, basic distribution of temperature and salinity is 
linear along the vertical and is not depend on time. The 
variables 6(t, x, z) and £(t, x, z) describe variations in the 
temperature and salinity about this main distribution. 
There are two types of thermohaline convection: the fin- 
gering ,7|, in which the warmer and more saline liquid 
is at the upper boundary of the area, and the diffusive 
type, in which the temperature and salinity are greater 
at the lower boundary In this paper we study the 
later case. 

The governing equations in the Boussinesq approxi- 
mation in dimensionless form are a system of nonlinear 
equations in first order partial derivatives with respect to 
time, that depend on four parameters: the Prandtl num- 
ber cr (usual value - 7.0), the Lewis number t (0 < r < 1, 
usually 0.01 — 0.1), and the temperature Rt and salinity 
Rs Rayleigh numbers 0, ^| : 



(d t - ctA)AV> + <r(R s d x Z - R T dJ 
(d 4 -rA)£- 5^ = •/(£,# 



= J(A^V), 



(1) 



Here the Jacobian J(f,g) = d x fd z g — d x gd z f is intro- 
duced. First equation in describes liquid particle 
pulse evolution in terms of stream function, second and 
third ones describe temperature and salt diffusion respec- 
tively. The boundary conditions for the dependent vari- 
ables are chosen to be zero, which implies that the tem- 
perature and salinity at the boundaries of the area are 
constants, the vorticity vanishes at the boundaries, and 
the boundaries are impermeable: 



ip = d*i) = 6 = £ = 0onz = 0, 1. 



(2) 



'Electronic address: skozi@poi.dvo.ru 



These boundary conditions are usually called free-slip 
conditions because the horizontal velocity component at 
the boundary does not vanish. 
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As a space scale the thickness of the liquid layer h is 
used. As a time scale value to = h 2 /x is used, where 
X is the thermal diffusivity of the liquid. Velocity field 
components are determined as v z — (x/h)d x ij) and v x = 
— {x/h,)d z ip. For temperature T and salinity S we have 
relations: 

T{t, x, z) = T_ + ST [1 - z + 6(t, x, z)) , 
S(t, x, z) = S- + SS [1 — z + £(t, x, z)] . 



T + - T_, SS = 



SY — S-, where T + , T_ 



Here ST = 

and 5+, 5_ are the temperatures and salinities on the 
lower and upper boundaries of the region, respectively. 
The temperature and salinity Rayleigh numbers can be 
expressed as follows: 



Rs = 



91 ■ 



-SS, 



XV xv 

where g is the acceleration of gravity, v is the viscosity 
of the liquid, a' and 7' are the temperature and salinity 
coefficients of volume expansions. 



we get a system of algebraic equations with solvability 
condition, having a form of a third order polynomial with 
respect to A: 



(A + ae 2 x 2 )(\ + e 2 x 2 ){\ + ts 2 x 2 ) 
+ — (\ + 1 e 2 x 2 ) =0. 



(5) 



Here x 2 — k 2 + n 2 n 2 is a full wave number and 7 is a 
constant: 7 = r + (1 — t)/N 2 . Equation JSJ) is known as 
dispersive relation and has three roots, two of which can 
be complex conjugates for a sufficiently small value of e. 
In the later case a Hopf bifurcation take place when at 
some values of N and e real part of the complex conju- 
gates roots turns to zero. It is true when 



e 4 < 



N: 



k 2 1 


'1-T N 


T 2 X & { 


1 + a , 


1 - T 


-e 4 — 
k 2 


l + (7 


k 2 (1 




X 2 ll 


+ a J 



e 4 T 2 x 4 . 



II. FORM OF CONVECTIVE CELLS AT LARGE 
RAYLEIGH NUMBERS 

Consider thermohaline convection in a limit of large 
Rs, which is true for the most of oceanographically im- 
portant applications (Rs w 10 9 — 10 12 ). After rescaling 
of the time t — (cri?s) _1 / 2 <', and the stream function 
ip = (aRs) 1 ' 2 ^' , we can rewrite basic system Q in a 
singularly disturbed form (primes are omitted): 

(d t - ae 2 A)A^ + {d x £ - (1 - N 2 )d x 8) = J(A^, V), 
(d t -e 2 A)0-d x ijj = J(0,iP), (3) 

Here a small parameter e 4 = 1/ <jRs and a buoyancy fre- 
quency TV 2 = 1 — Rt/Rs are introduced. In this system 
singular perturbations are present as e 2 before Lapla- 
cians. If we let e — 0, then our system © turns into 
common equations, describing two-dimensional internal 
waves with the constant buoyancy frequency N in the 
Boussinesq approximation. 

For investigating of a linear stability problem for the 
system with boundary conditions (J2J we omit nonlin- 
ear terms in the right part of the system and choose a 
solution in a form of normal mode: 



ip(x, 


z,t) 


= Ae xt ~ 


lkx sin rnrz, 




0(x, 


z,t) 


= a T e xt 


~ lkx sin nnz, 


(4) 




z,t) 




~ lkx sin nirz. 





where A is an eigen value, describing growth rate of the 
mode, k is a horizontal wave number, n is a number of 
the mode and A is an amplitude of the mode. After 
substitution of the expressions (0J into the system © 



where w = Im(A) is a Hopf frequency. 

Because dispersive relation (|SJ) explicitly contains small 
parameter, we can choose one of the complex conjugates 
roots and express A in the form of an asymptotic expan- 
sion by the powers of e: 

A = A + e 2 Ai + e 4 \ 2 + e 6 X 3 H 

After substitution of this expression in (jSJ) we have for 
A,-: 



A a — 

A 2 = 

Fi = 

F 2 - 

F 3 = 



-(k 2 /x 2 )N 2 
-(x 4 /X )F 2 , 
(7-d)/2>0 
(3F 2 + 2dFi - 
4F? + 4&F? + (Cl + C 2 )F 1 
+ (t + <j)(C 1 +t<t)/2 > 0, 



X 1 =x 2 F 1 , 
A 3 = -x 8 /(k 2 N 2 )F 3 , 

C 2 )/2 > 0, 



where constants are: C\ = 1 + r + a and C 2 = r + 
fi + rcr. It is interesting to note, that when N* > N > 
functions F±, F 2 , F 3 are positive. The growth rate, caused 
by thermohaline convective instability can be written as 
follows H3: 



Re (A) = e 2 x 2 F 1 - e e 



k 2 N 2 



F, 



(0) 



One can see (FIG.^), that for a given mode with number 
n, the growth rate is maximal for some k, which deter- 
mines horizontal size of the most prominent convective 
cells. Also the first convective mode has maximal growth 
rate, so that convective cells tend to be tall and thin. 
Really about 30 first modes for N — 0.3 have positive 
growth rate (FIG. |2J, but the most prominent cells any 
way will be tall and thin, as having relatively more large 
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FIG. 1: Grows rate Re(A(fc)) (a), and frequency Im(A(fc)) (b) 
of travelling waves Q for the first convective mode. Here are: 
e = 0.00153, a = 7,t = 1/81 and N = 0.3. Curves 1 are one 
term approximations; curves 2 are two term approximations; 
curves 3 are exact solutions of the equation J5J . 
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FIG. 2: Grows rate Re(A(fc, n)) for convective modes with 
different mode number n. All parameters are the same as in 
the FIG. [H 



rate of growth. For simplicity further we consider only 
the first convective mode. 

Although above developed perturbation approach 
gives qualitatively true estimates, for more accurate re- 
sults one should immediately solve algebraic equation JSJ. 
Rewrite it in another form, with new introduced variables 
P 2 = e 2 x 2 « e 2 k 2 , X = \/P 2 and Y = N 2 /P 4 . 

{X + <j)(X + 1)(X + t) + Y(X + 7 ) = 0. 

Roots of this equation depend on parameter Y", so that 
finally A depends on horizontal wave number P. 

Consider actual oceanographical system such as an in- 
version of thermohaline staircase. Let it has thickness 
h = 250 cm and temperature difference ST = 0.1°C, 
also (j — 7 and r = 1/81. In this case e = 1.53 x 
1CP 3 and non-dimensional critical buoyancy frequency 
AT* = 0.35136. For N = 0.2764 the most unstable 
mode has P* = 0.1599579 and width of convective cell 
l c = neli/P* « 7.7 cm. For comparison formula © gives 
P* = 0.126, i.e. somewhat less than the exact value. 
From @ one can extract dependence of P* from N, hav- 
ing form P* = [A^Fi/^Fs)] 174 - 

From picture (FIG. |SJ one can easily see, that when 
value of N becomes slightly small than its critical value 
AT* w 0.35136, value of P* abruptly (as (AT* - N) 1 / 4 ) 
increases and becomes maximal at A^ ss 0.3. When A^ 
becomes even more less P* decreases to P* w 0.13688 
for N = 0. It should be emphasized that P* is nearly 
independent of A^ in considered case. 

This result for thermohaline convection at large Rs is 
sufficiently different from that for small Rs, when criti- 
cal wave number is K=n/V2 [13. For our case typical 
wave numbers are of the order 0.1/e. This estimate is 
more accurate, than mentioned in pj. Thus convective 
cells for large Rs have tall and thin geometry from linear 
stability analysis. Physically one can understand this ef- 
fect from consideration that when Rayleigh numbers are 
large, the buoyancy forces acting vertically are also large 
in comparison with forces of inertia of liquid particles 
determining width of cells. 

III. FINE STRUCTURE GENERATION AND 
DIFFUSIVE CHAOS 

Now we will study nonlinear vertical modulations of 
amplitude of our tall thin convective cells. 

At first, introduce new small parameter, extracted 
from geometry of the convective cells: e = l c /h — 
ne/P* » 20e. Rescale variables tp = e 2 tp' 7 9 = e6' , 
£ = e£' (prime will be omitted), and introduce one more 
small parameter E = e / e ~ 1/20. After changing of the 
space scale from h to l c basic system @ appears in the 
form: 

(d t - aP 2 A)A^ + (d x £ - (1 - N 2 )d x 6) = J(A</>, V), 
(d t -E 2 A)8-d x iP = J(0,il>), (7) 
(d t -TE 2 A)£-d x ^ = JfoVO. 
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FIG. 3: Dependence of the wave number P» = ke of the most 
unstable mode from buoyancy frequency TV. Curve 1 is the 
exact numerical solution; curve 2 is the approximation by 
formula JSJ. 



At second, introduce slow vertical variable Z = ez and 
slow time T — e 2 t. In accordance with multi-scale 
decomposition technique we get prolonged derivatives 



aecompositi 
dill El: 



d z ed Zl 

d t -> d t + e 2 d T , 

A - 9 2 + e 2 t?i, 

A 2 - ^ + 2 e 2 9 2 a 2 + e 4 9| ; 

ftA -> 9 t a 2 + e 2 9 T <9 2 + e 2 a t 5|. 

Let buoyancy frequency somewhat less than its critical 
value N 2 = N 2 — e 2 R. Parameter R here is a forcing of 
the system. Equations Q now get a form: 

& - <rE 2 d 2 )d 2 ^ + (dxt - (1 - N 2 )d x 9) = 
~eJ z ^,d 2 x ^)-e 2 [(d T d 2 x +d t d 2 z 
-2aE 2 d 2 x d 2 z )^-Rd x 9], 

(d t - E 2 d 2 x )e - d x i> = 

-eJz{i>i6)-e 2 {d T -E 2 d 2 )6, (8) 

{d t - rE 2 d 2 x )i - d x1 p = 

-eJzty,® - e 2 (d T -rE 2 d%)i. 

Solutions of these equations we will find as the asymp- 
totic sets by powers of the small parameter e: 



ip — eipi + e ip2 + e ip 3 + 
9 = e6 t + e 2 6 2 + e 3 9 3 + ■ ■ 
£ = e& + e 2 & + e 3 & + • • 



(9) 



After substitution of these expressions into equations (JSJ 
collect terms at the same powers of e. As a result we 
have systems of equations for determining of the terms 
of the sets ©. Thus, at e 1 we have following system: 

{dt - <jE 2 d 2 x )d 2 ^ + (cUi - (i - Nl)d x e x ) = o, 

(d t - E 2 d 2 x )6 l - d x i>t = 0, (10) 



Choose for this system solution in the form of normal 
convective mode travelling to the right, with constants 
of integration B(T,Z) and C(T,Z), depending on slow 
variables. 



ipi = A(T,Zy wt - lKx + c.c. 

0i = a T (T, z)e iuJt - iKx + B(T, Z) + c.c. (11) 

£i = a s (T, Z)e iut - iKx + C(T,Z) + c.c. 

Here wave number K = (N) is a horizontal wave num- 
ber, corresponding to the most unstable waves of convec- 
tion, and maximal value of K = n from choice of the 
space scale, related with convective cells. It is attained 
when N ss 0.3. Parameters of the normal mode lllfl are 
related as follows: 



iK 



Qt = 



iK 



icu + E 2 K 



as 



tE 2 K 2 



A, 



(iuj + aE l K 2 )(iuo + E 2 K 2 )(iuj + tE 2 K 2 ) 

+N 2 (iu + 1 E 2 K 2 ) = 0. 

Last formula is actually the dispersive relation (J5J, but 
for an infinitesimal vertical wave number. Also for critical 
buoyancy frequency N*. and wave frequency ui we have: 



n: = 



1 - T 
1 + cr 
1 — T 



4 7^4 



- (1 +t){t + a)E^K 
t 2 E a K\ 



1 + cr 

oj 2 = N 2 + (a + T + ar)E 4 K 4 . 

System of equations at e 2 is the same as (|10|l and does 
not lead to any new results. System at e 3 is: 

(d t - crE 2 d 2 x )d 2 x ^ 

+(dM - (i - N 2 )d x e 3 ) = -Jz^d^i) 

-(d T d 2 x + d t d\ - 2aE 2 d 2 x d 2 z )4>i + Rd x 6 x , 
(d t - E 2 d 2 x )9 3 

-d x ^ 3 = -J z (i> u ei) - (dr - E 2 d 2 z )6 1 , 
(d t - rE 2 d 2 x )i 3 

-d x ^ 3 = -Jz(ik^i) ~ (dr - TE 2 d 2 z)ix. 

After substitution into the right parts of these equations 
expressions <|11|) we get a system with resonating right 
parts, breaking regularity of the asymptotic expansions 
© ■ The condition of the absence of secular terms in this 
case takes form of so called ABC system |(j (intermediate 
calculations are omitted): 

d T A = E 2 /3 1 d 2 z A + R/3 2 A- /3 3 Ad z B + /3 4 Ad z C, 
d T B = E 2 d 2 z B-E 2 Mz\A\ 2 , (12) 
d T C = TE 2 d\C -TE 2 [3 6 d z \A\ 2 . 

Here coefficients are: 
1 



(dt - rE 2 d 2 x )Ci - c^Vi 



0. 



1% 



1 



[(iuj + oE 2 K 2 ^ 



E 2 K 2 



(1 - r)E 2 K 2 
(ico + rE 2 K 2 ) 2 



( iu 



2cr 



\E 2 K 2 

(iuj + gE z K 



iuj + E 2 K 2 

, (l-r)iw 



(ioj + tE 2 K 2 ) 2 



fa 

05 



{iw + oE 2 K 2 ) + 
Po 1 



1 



iu + E 2 K 2 
27y 4 



iu + tE 2 K 2 
04 = 



Po 



0o 



E 4 R 4 ■ 
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iuj + rE 2 K 2 ' 
2K 4 



l E A K A 



Thus in this article we have derived ABC system of am- 
plitude equations for travelling waves of double-diffusive 
convection in a limit of high Hopf frequency (large R$) 
in the infinite horizontal layer. 

Equations i|12|) have nontrivial solutions, describing 
such phenomena as diffusive chaos and can be used 
for simulation of formation of patterns, like vertical fine 
structure of temperature and salinity in some areas of the 
ocean, for instance, in inversions of thermohaline stair- 
cases. 

Transform system (|12[) to more convenient form by in- 
troducing new time variable T" = E 2 T, and applying 
following substitutions: 



A' = AE^y/fclfclexpi-iRfcRT'E- 
B' = \p 3 \E- 2 B C = \(3i\E- 2 C. 



System <|1 2|) now gets form (primes are omitted): 

d T A = (3 1 d 2 z A + a 2 RA- a 3 Ad z B + a i Ad z C 1 
T B = d 2 z B-d z \A\ 2 , 
d T C = rd 2 z C -Ta 6 d z \A\ 2 . 



(13) 



Here are coefficients: a 2 = P 2 rE 2 , «3 = 03/\03\, a 4 = 
4 /\0 4 \ and a 6 =0 6 \0 4 \/(0 5 \0 3 \). 

We developed numerical models for parallel calcula- 
tions of system Ijl3|) based on explicit and Dufort-Frankel 
schemes. For numerical experiments it were chosen pa- 
rameters a = 7, K — 7r, E — 0.05, and two values of 
Lewis number: r = 1/10 and r — 1/81. Governing pa- 
rameter R was in the range from 0.1 till 50. Number 
of vertical dots n was from 256 till 2048 to resolve mi- 
crostructure. 

System (I13[) was integrated numerically on multipro- 
cessor computer MVS-1000/16 with zero boundary con- 
ditions and sinusoidal initial conditions for dependent 
variables. In the most cases the initial state was de- 
stroyed after some time via a multiple Eckhaus instabil- 
ity (birth of convective cells, 01) an d was followed by 
diffusive chaos state, with strong space-time irregularity. 
In this case mean profiles of temperature and salinity 
become perturbed so that all layer of inversion becomes 
divided on 10 - 30 small layers (see FIG. H for r = 1/81 
and time t = 9.28 hours). Buoyancy frequency (FIG.JSJ) 
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FIG. 4: Numerical solution of the system (1131 at the time 
t = 9.28 hours. Here variables: \A(T,Z)\ (a), B(T,Z) 
(b), C(T,Z) (c). Forcing parameter is R = 16, also e = 
0.00153,a = 7,t = 1/81. Layer depth h = 250 cm, the 
number of grid points is 2048. Dimensional variations in tem- 
perature and salinity are proportional to B and C, its can be 
estimated as 6.9 • 10~ 5 ■ B ■ ST°C and 6.0 • 10~ 5 • G ■ 5S% 
respectively. Amplitude of stream function is proportional to 



|y4| and can be estimated as 6.9 • 10 
variation in temperature, for instance, 



\A | cm /sec. Maximal 
is about 20% from ST. 
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FIG. 5: Buoyancy frequency N (cycles per hour) vertical 
microstructure for time t = 9.28 (a), and t — 11.6 (b). Other 
parameters are the same as in the FIG. HI 

becomes very irregular, and all this fine structure slowly 
changes with the time. 

IV. CONCLUSION 

In this article we developed mathematical model, de- 
scribing formation of vertical convective patterns in two- 
dimensional double-diffusive convection in a limit of high 
Hopf frequency for an infinite horizontal layer. A phys- 
ical system, corresponding to such model is inversion of 
thcrmohaline staircase. Some typical parameters of in- 
versions are presented in the tabled It is known that 
parameters of stratification in the inversions are often 
near the onset of convection. Also vertical microstruc- 
ture (usually step-like) are often observed in the inver- 
sions along with small scale turbulence [lj(. Results of 
this work are in qualitative agreement with these obser- 
vations. For more comparison see also [j| and references 
therein. 

Although we aimed to construct mathematical model 
with hne structure generation in double-diffusive system 
without its detailed relation with experimental data, we 
should note a few points of such relation: 



TABLE I: Parameters estimations for inversion of thermoha- 
line staircase. For all cases it is true T = 15° C, S = 36% , 
a = 7, t = 1/81, to - diffusive time scale, t' - main time scale, 
Nq = (1 — t)/(1 + a) - limit of critical buoyancy frequency 
at e = 0, 400 to ~~ time of establishing of diffusive chaos in the 
inversion. 



Parameter 


1 


2 


3 


h [cm] 


400.0 


250.0 


100.0 


ST[ C] 


1.0 


0.1 


0.1 


5S [% ] 


0.33 


0.033 


0.033 


Rt 


9.5 x 10 11 


2.3 x 10 10 


1.48 x 10 9 


Rs 


1.08 x 10 12 


2.6 x 10 10 


1.69 x 10 9 


e 


0.0006 


0.0015 


0.003 


e 


0.012 


0.03 


0.06 


l c [cm] 


4.7 


7.7 


6.0 


to [sec] 


1.12 x 10 s 


4.38 x 10 7 


7.0 x 10 6 


t' [min] 


0.68 


1.7 


1.07 


N [eye/ hr] 


4.95 


1.97 


3.13 


400 4 [hr] 


4.5 


11.4 


7.14 



• Developed model predicts that fine structure 
should exist in given system for a wide range of 
parameters with typical time of pattern formation 
of about a few hours. 

• Fine structure has very irregular shape, slowly 
changing with the time in accordance with solution 
of ABC-system of diffusive chaos type. Mean pro- 
files of temperature and salinity become perturbed 
so that all layer of inversion becomes divided on 
10 — 30 small layers. 

• Buoyancy frequency structures have peak emissions 
with amplitude and width with reasonable agree- 
ment with observational data 0] in the cases when 
temperature and salinity differences per layer are 
relatively small, as in the table [I] 

Of course, developed model, based on ABC-system of 
amplitude equations has its limitation of weakly non- 
linear approximation, which hardly allows to get "full- 
fledged" step-like vertical structure of density, as noted 
in p • Partially this defect can be overcome by regarding 
amplitude equations arising at higher orders of small pa- 
rameter in multi-scale decomposition method. But this 
is a matter of further works along with more detailed 
comparison of predicted fine structure parameters with 
experimental and observational data. 
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